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Ce eR 1) 


INTRODUCTION 


Pel Our Changing Environment 


More technological advances have been made in the Twentieth 
Century than in all of man's past history. These advances have affected 
nearly every aspect of existence on this planet as well as fulfilling the 
age-old dream of setting foot on other heavenly bodies. Man, as a result, 
enjoys a longer and easier life. 

These phenomenal advances have not been made without some 
potentially harmful side effects to the environment. Factories have con- 
tributed to air and water pollution. Questions are being raised about the 
relative merits of Supersonic travel versus the adverse effects of the 
resultant sonic boom. Every man has a car inhis garage but the country- 
side is being ruined by ribbons of concrete. Many recent discussions have 
centered on the adverse effects of electromagnetic and particulate radiation 
on biological systems in general and human tissue in particular, but there 
are also many advantageous uses for electromagnetic and particulate 
radiation such as the use of electromagnetics in diathermy treatment of 
many musculoskeletal disorders and the use of nuclear energy to generate 


erectimlcret@e l/c, 


1.2 Heat for Therapy 


Man has searched for more effective methods of applying heat to 


his body for various aches and pains since he first discovered what the 





heat of the sun could do for him. This search has led man to investigate 
nearly every method of heat production including the sun, hot water, 
infrared radiation, ultrasonics and electromagnetic radiation. In recent 
years, there have been many meee’ of successful applications of electro- 
magnetic heating of deep body tissue for therapeutic purposes (diathermy). 

Since diathermy has proven to be effective for treatment of many 
musculoskeltetal disorders, the search is continuing in the Department 
of Physical Medicine and Rehabilitation at the University of entre core oa 
in an effort to find the most efficient method of realizing deep tissue heating 
without damaging the skin and Subcutaneous fat. 

Computer solutions of the mathematical models of the diathermy 
problem have in the past been confined to investigating only a rather 
limited set of points in the volume of the models. This research will 
investigate an algorithm which can more efficiently provide a much more 
complete picture of the heating. 

A direct consequence of diathermy research will be the develop- 
ment of the analytic methods which can be uSed in the study of the adverse 
effects of electromagnetic radiation and as a quantitative aid in the estab- 


lishment of radiation safety standards. 


Daomewadiation Safety Standards 


The relative importance of determining the biological effects of 
electromagnetic radiation can be ascertained by the attention given to 


the problem in recent years: 





a. Numerous congressional hearings were held on the subject in 
1967 and 1968. oa 

b. The Radiation Control for Health and Safety Act of 1968° 
(RCHSA-1968) was enacted as a result of the congressional 
hearings. 

c. The Armed Forces have, since 1957, held an annual Tri- 
Service Conference on the Biological effects of Microwave 
Radiating Equipments, Z 

d. There is continuing research into the possible uses of electro- 
magnetic radiation in the field of medicine. 

e. There has been a great deal of research conducted on the use 
of electromagnetic heating of deep body tissue for therapeutic 

purposes (diathermy). capt 0 

Discussion and research into the biological effects of electromagnetic 


radiation have resulted in both the United States and the USSR accepting 


radiation safety standards for prolonged exposure to radiation. 


1.4 The Need for Research 

That there is a need for further research can be evidenced by the 
fact that the safety standards set by the two countries differ by three orders 
of magnitude. In fact, the RCHSA-1968° directs the Secretary of the 
Department of Health, Education and Welfare to conduct studies of the 
problem and report the findings to Congress. 

How should the research be conducted? Inthe past, most of the 


research on the biological effects of electromagnetic radiation has been 





qualitative rather than quantitative. Considerable research has been 
conducted on animals and the results have been extrapolated to man with 
a fair amount of success. But infact, the very question of qualitative 
versus quantitative research has been debated at great lengths while the 
proponents of each side continued to conduct the research as best they 


could in view of the facilities at their disposal. 


1.5 An Engineering Approach 


A reasonable engineering approach to the problem would incorporate 
both quantitative and qualitative research in a team effort by engineers and 
physiologists where members of each discipline would investigate that 
aeieect ot the problem to which they are best suited. 

The engineers are well suited to investigate the quantitative aspects 
of the problem by first constructing planar and cylindrical layered models 
which exhibit the electromagnetic properties of human tissue. This analy- 
sis would mainly encompass the problem of determining by analytical 
methods the electric and magnetic fields within the tissue layers of the 
models for different internal or external source configurations. 

When the models have been developed to the point where they reliably 
agree with theoretical calculations, the same modelling materials can be 
used to model the complex geometries of the human body. The answers 
derived from the models can then be correlated with data from qualitative 


methods. 





1.6 Research Objectives 

The approach outlined in the preceding section is currently in use 
at the University of Washington. The problem under investigation is the 
determination of the optimum frequency of electromagnetic radiation and 
applicator configuration to be used for therapeutic heating of deep tissues. 

The subject of this particular research is a comparison of the 
Guassian quadrature and fast Fourier transform (FFT) algorithms used 
to compute the electric fields induced in planar layers of human tissue. 
The vehicle problem is the computation of the relative heating pattern 
arising from a direct contact diathermy applicator. The applicator has 
a TE V0 waveguide mode radiating into a bi-layered, semi-infinite lossy 
medium (See Figure 1.1). A subcutaneous fat layer overlying a semi- 
infinite muscle is modelled by this medium. 

Current analysis utilizes the Gaussian quadrature method to per- 
form the necessary transformations. The economics of computation 
time in the continuing research program dictate that only a relatively 
few points in any plane of the tissue can be computed and those must be 
along an axis of symmetry. The fast Fourier transform will compute 
an entire array of points in any desired plane, thereby presenting a more 


complete picture, at a considerable savings ona cost per point basis. 
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1.7 Content of Succeeding Chapters 


In Chapter 2 the electric field equations and relative heating 
expressions to be solved are developed from the literature. Chapter 3 
discusses the fast Fourier transform which will be used to accomplish 
iemiveroctranorOrms required to compute the electric fields for the 
relative heating patterns. 

The use of the FFT in this and another related problem will be 
presented by discussion and flow diagrams in Chapter 4. 

This research project analyzes the merits of using the FFT to 


compute the electric fields necessary to obtain the relative heating patterns. 








The results from the fast Fourier transformations at several frequencies 
are compared in Chapter 5 to results obtained from Gaussian quadrature 
transformations already accomplished by Guy and McDougall. 

A brief discussion of the conclusions to be drawn from this 


research is presented in Chapter 6. 





Crago Re 2 
Pees tiVE HEATING IN’ PLANAR LAYERS OF HUMAN TISSUE 


pee Objectives 

Relative heating from a rectangular aperture in an infinite conducting 
plane radiating into a two layered model of human tissue is defined in Sec- 
fton 2.2 (See Figure 2.1). The human tissues modelled are a subcutaneous 
fat layer of depth zy Overly Meet Seinil—iniimite muscle layer. The succeed- 
ing sections develop the necessary electric field equations required for 
computation of relative heating. 

The development of the field expressions follows the approach of 
falienicuver - and Guy. = Villeneuve utilizes the Fourier transform 
approach in computing the effective input impedance of an open-ended 
waveguide radiating into a plasma environment. This problem arises in 
designing communication antennas on high speed reentry vehicles used in 
space exploration. 

Villenetime wim his work, develops general expressions for the fields 
exterior to the waveguide; then, he looks at the expressions only at the 
aperture. Guy extends the expressions of Villeneuve by matching boundary 
conditions and derives the field expressions necessary to compute relative 


heating. 


2.2 Relative Heating 


The criteria for the design of an optimum diathermy applicator 


for treatment of musculoskeletal disorders is that the heating per unit 
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FIGURE 2. 1 


Tissue Model and Coordinate Definition 
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volume in the subcutaneous fat be significantly lower than the heating per 


unit volume in the musculature. 
Oicecriteria. Relative heating at any point (x, y, z) 
eeune per unit volume at the point (x,y, z) 


unit volume at the point (0,0, Z,), where z 


3 
cutaneous fat layer. Guy 


Boespoint (x,y,z) is given by 


WT 


] 


Ww +w + W 
x y Z 

0) 

o ie f 

5 

—|E 

5 lef 


SE f 


The term relative heating arises from 

is defined as the 

Peltiihve tourmle Neatine sped 
is the thickness of the sub- 


has shown that the heating per unit volume at 


where 0 is the conductivity of the medium and the E's are the electric 


field strengths at the point (x,y,z). 


Medting at a point (x, y, 2) 
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=e 


aicrintae muscle by 


Pp i 
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in the fat layer is given by 
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From these expressions the relative 
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e+ ef +e 


Pp = Leo 
a f Z 2 

m Eo mee +|E | 

where 
= E (0 ao 

E 0 Bi , 0,2)) 232 69 
E . = E (0,0,z.) 2a 
y9 y l 
Eo = E (0, 0,2,). Ziel 


The W (0,0,z ) which appears inthe p and p expressions 

a 1 nt Ay 

m f m 

fenmiaiizes the relative heating to unity at the fat-muscle interface. 
Ewoepaeneral Solution of Electric Fields in Layered Media 

In rectangular coordinates, general solutions of the scalar wave 

. : i 10 

equation may be obtained as a linear superposition of plane wave solutions. 
This process leads to the development of a complete set of Fourier trans- 
form pairs which completely specify all of the field components throughout 
the volume of the fat and muscle. 

The electromagnetic fields may be derived from two parallel 
vector functions whose direction in space is fixed. Assume these func- 
tions to be 
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e m 





eZ 


a + ee w= Q. Ze 


A linear superposition of plane wave solutions apropos to this 


Peoolerm (see Figure 2.1) is, for Region (1) 
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From Maxwell's equations 
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iS 
Making the appropriate substitutions and carrying out the indicated 
operations yields the following equations for the electric fields (See 


igure 2.1). Region (1) 
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Matching boundary conditions at z= 0 and z= zy 


and solving 


the six resulting simultaneous equations gives the various f's and g's. 


This has been done by ony with the following results for the fat-muscle 


case under investigation. 


Region (1), in the fat 


CO 
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\ Voc vez) Bee : ae ay 
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00 
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mopton (2), in the muscle 
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= Zz 2 2 re 
Bt = (w+ v°)(W 2, (ea Meaewak « oes 
l Z 4 l 
4: 2% 2.2 % 
= k se ; oir 
C aS €1 Ww, k 1 2.3.34 


Since the source field is at z= 0, then forall z> 0 the fields 
must represent waves travelling and attenuating in the z-direction. That 
is, the branches of the propagation constants Wi and Ww, must be 


chosen to yield these characteristics. Such a choice is 








Zz s 
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Za 2 a4 Z 3! Z 
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Now, the & which appears in each of the electric field equations 


is the Fourier transform of the electric fields at aperture, l.e., 2 = 0. 


2.4 Electric Fields in Single Semi-infinite Layer 
When only a single semi-infinite muscle layer is considered, the 


equations in the previous section can be completely utilized by setting 
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z,. equal to some convenient value such as one (1) centimeter for 


i) 
4 6 
a6 “1 


= VW 6. Se 


normalization purposes and by letting Wi > 1 2° 


Another method of obtaining the expressions for a single layer 
medium is to utilize the equations for region (2) and match boundary 
conditions at z= 0. Both methods are easy to carry out, and since the 


results are relatively simple, they will not be included in this discussion. 


2.5 The Theoretical Problem 
An independent choice of any two of the field quantities will auto- 
. Te EZ 
matically determine all other field quantities. For purposes of this 
research, the field quantities will be chosen to describe a field configura- 


tion at the aperture which is that corresponding toa TE waveguide 


10, 


‘ 


AOC . 

The electric field in the x-direction is chosen to be of a uniform 
magnitude and phase along the x-direction and a cosine magnitude varia- 
tion with constant phase in the y-direction. The electric field in the 
weairection is zero everywhere in the aperture. 


This can be mathematically stated as 


b 
= -x< wae aes 2 
E(x, y, 0) 0 oc<cy < : 5 y Be) 
E(x, veeuy = Cos = a 5 5 
coum) = © -O< XS © Ze reiee 


18 
2.6 Errors Due to Less Than Infinite Integration 

A completely accurate computation of the electric fields at any point 
(x,y,z) would necessitate integration over the complete infinite u- v 
Plane. Such an analysis in this problem is entirely unrealistic and would 
usually make only extremely insignificant contributions outside the normal, 
finite ranges of integration. 

Occasionally, for some integration methods, significant computa- 
tional errors will result when the integration is carried out over finite 
limits. For the problem defined in the previous section, two Specific 
instances exist for which finite integration fails to yield answers which 
are devoid of significant computational errors. 

First, when electric fields exist near a conducting wedge (See 
Figure 2.2), a phenomenon known as edge effect arises. Near a con- 
duction wedge, some field quantities will tend toward infinity as the edge 
is approached even though the energy remains finite. This edge effect 
will produce very high field strengths near the edges of the aperture and 
hot spots in the relative heating pattern. 

In the analysis of the edge effect, Gun assumed a single homo- 
geneous semi-infinite lossy medium. This assumption does not account 
for reflections from the fat-muscle interface, and when Zs the distance 
from the surface to the fat-muscle interface, is small, this assumption 
introduces another error which must be corrected. 

ee has analyzed the two problem situations described above and 


derived corrections which are added to FE. and ET when x and z are 
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iieclose proximity to the aperture edge or when zy iS so small that 


imaging effects must be taken into account. The corrections are given by 
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These corrections are to be added to the electric fields at points 
near the edge of the aperture and at all points in the Region (1) expres- 
sions for E.. and eS when zy is small enough to cause imaging. Note 
Peetaeeinteprands in the correction can be subtracted from the integrands 
Semeeceim ocction 2.3 before the integration with respect to u is performed. 

All of the necessary equations for computing relative heating pat- 
terns have now been developed. Attention can now be focused on an 


efficient algorithm to compute the necessary transformations. 
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THE FAST FOURIER TRANSFORM 


B. 1 Background 


Since 1822 when Jean Baptiste Fourier (1768-1830) published his 
imeorlve Analytique de la Chaleur, in which he set forth the trigonometric 
series and integrals which are named after him, Fourier analysis has 
played an important role in the analysis of many of the problems of physics 
and engineering. Fourier himself only presented the general boundary 
value problem for the case of heat conduction, but today his methods find 
wide application in the disciplines metnioned above. Besides the heat 
conduction problem where the method is still apropos, many problems of 


today are particularly amendable to the method; to name a few: 


a. Power spectral density 

D: System simulation 

GC. Noise analysis 

d. Statistical distribution nea ove 

ee Frequency analysis of bandlimited waveforms 
i Convolution (lagged product) problems 

g. Acoustic wave propogation 

he Speech transmission 

ie Transport phenomena 

ie Optics 


k. Mlectro magnetics. 





Ne 
3.2 Fourier Analysis 
The theory of Fourier analysis is beyond the scope of this 
£215 ied . 
report, but several excellent texts exist pertaining to this subject. 
The major drawback of Fourier analysis lies in the need to evalu- 


ate integrals which have infinite limits, as is seen in the one-dimensional 


Fourier transform pair utilized in frequency analysis: 


CO 
S(f) = \ eae ian 3.2.1 
=- © 


_JeT be 


C—-> 8 


mit) = 5 s(£) df 3.2.2 


8 


or the equivalent two-dimensional form usually seen in boundary value 


problems: 
00 
l -j(ux + vy) 
Siu. v) = ES d(x, y) e dxdy bee: 
t - 0 
00 
] j(ux + vy) 
ioc) a ae S(u, v) e dudv 3.2.4 
= 00 


Equations 3.2.3 and 3.2.4 are typical of those developed in Chapter 2, 
but for purposes of illustration of the FFT the equations 3.2.1 and 


3.2.2 will be used. 


3-3 The Fast Fourier Transform 
Prior to the advent of high speed computers, complex integration 


was handled by such methods as the saddle point or stationary phase 





Zo 
techniques which are valid only for expressions of large argument. Thus 
the really important questions of analysis of near zone fields of an antenna 
were seldom investigated since in this area the integrals were unassailable 
by analytic methods, except ina few special cases. The computer changed 
this, but did not significantly advance the state of the art due to the high 
cost of computing time. 

A real advance came when Cooley and Tukey! found the algorithm 
known as the fast Fourier transform. It then became economical to com-- 
pute discrete Fourier transforms, and many problems could be economi- 
cally investigated. 

Returning now to the Fourier transform itself, and in particular 
looking at the discrete Fourier transform pair associated with equations 


meer leeand 3.2.2, 


Nat -jon f te N 
Sf.) = At ) x(t) « Mea Oct Lee Oral 
n k 2 
k=0 
N/2 jen fk 
x(t, ) uN E » S(f je he= O74 tZ,2...N-1 3x32 
n=N/2 
now, making the following substitutions, 
tL = kAt t = T/N 
f =ndAf i= es 
ie) ak 


then equations 3.3.1 and 3.3.2 become, according to Cooley, 
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N- 1 
-j2 k)/N Aa 
S(n) =) x(k) © J" (KI/ > 
k=0 te-sO,b ete, 4 .. N=] 


Ole. Nal 
N-1 
x(k) = P S(n) elém (nk)/N anon 
bea 


3.4 Factoring the Transform 


For purposes of illustration, when equation 3.3.3 is written out 


Lo Jer k/N 


in full for N=4 andwith W Brigham and videronae 





define 
w° w° w° W x (0) S(0) 
w° w w° W ae) S(1) 
[S(n)] = ee = S47 al 
w WwW WwW w x4 (2) S(2) 
w? w? w° w x 9(3) S(3) 


——— SS — = 


Now define [S(n)] as: 





S(0) X (9) 
S(2) x 9(1) 

[ S(n)] = = = [W] [xp (n)] 
S(1) Xo (2) Cy 
S(3) x 9(3) 


0 
The [W] matrix in 3.4.2 can now be factored, and since W = 1 





ZS 
0 0 





Pew OO ii w § «O X (0) 
1 w 0 0 lao we X(t) 
o(m)} = 1 3 3A 
0 0 1 W ieOa vy (0 X (2) 
00 1w} foro wy |x (3) 
Alternately, this can be written: 
1 we 0 0 
1 w- 0 0 
metn)) = [x (n)] = 1 [x,(n)] Sr 
2 00 1wW 
0 0 I w? 
where: 
1 0 ae 0 
Om lad we 
ieca(n)| = [ x(n) ] 3.4.5 
, iG w* 0 : 
Om leo we 


This factorization looks somewhat mystical, so to put minds at 
ease and insure that the factorization is valid, and to insure that equation 
feomrteagrees with 3.3.3, from equation 3. 4. 4, f 

0 

x, (9) st x, (1) W 

2 

x, (0) 1F x, (LW 

etn) | = \ 3. 4.6 
Ne + (3) W 


x, (2) + x, (3)W ; 


and from equation 3.4. 5, 








ZG 








x (0) + 
Xo (1) + 
[x (n) ] = Dea e 
x (9) + 
X_(1) + 
then, substituting 3.4.7 into 3.4.6 yields 
% (0) + x 9(2)We + (x (1) + x9(3)W°)W" 
x (0) + x (2)w° + (x (1) + x yw yw? 
0 0 0 0 
[ S(n) ] = 2 2 1 3.74.16 
% (0) + X\(2)W + (x(1) + x \(3)W )W 
Z CaS 
x, (0) + (2) W + (x) (1) + X (3) W )W 
n+4 n+4k 


n 
If we note that W = W or, more generally speaking w= W 
Mieere kK iS Some integer, then 3.4.8 upon rearrangement becomes 


once again equation 3.4.2 written out in full, instead of matrix form: 


0 0 
X_(1)W + x) (2)W + x 


+ 


x 9(0) 


x 9(1)W + x 9(2)W t X (3) W 


+ 


x (0) 





+ 


x (1)We + x (2) W* + x_(3)W 


X (0) + x (1) W> + x 9(2)W° + X_(3)W 


x (0) 





One discrepancy has arisen in that we have derived [ S(n) ] instead 
of S(n). This can very easily be corrected by an algorithm which looks 
at n as a binary address of a particular element of the matrix S. To 
find the correct address is simply a matter of bit reversing the binary 


number n. 





oe 


Then, 
S(00) 
S(01) 
[S(n)] = 20 
S(10) 
S(11) 


and bit reversing the address yields 


S(00) 


S(10) 
[S(n)] = Bd ail 
S(01) 





S(11) 


While it is true that the above development applies only to the 
ease N=4 and not to the general case, when the discrete Fourier 
transform, equation 3.3.3, is examined for a general case, as has been 
shown by Gentleman and Sande, ae then the factorization can be performed 
for any N that is not a prime number. This factorization is a recursive 
Formula Which defines larger Fourier transforms in terms of smaller 
transforms. A multi-dimensional discrete Fourier transform, suchas 
this research investigates, can be efficiently computed by performing 


the factorization in each dimension separately. 


3.5 The Algorithm 


The obvious question to ask now, is how to obtain the factored 


iMeatmices Indicated in equation 3.4.3. This is most easily seen ina 
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Peee praph (See Figures 3.1 and 3.2) where a vertical column of nodes 
will correspond to the x) (Kk) arrays (where for our example ¢£ = 0, 1,2). 
For cases where N = 2%(y an integer), we willhave y+ 1 vertical arrays 
and N horizontal arrays. For this example we willuse y=2, N= 2° = 4. 
In the tree graph, a solid line will indicate that the node in the previous 
vertical array is to be multiplied by the W factor raised to the power 
indicated in the node of interest, then added to the value of the node from 
which the dashed line originates. The following algorithm yields the neces- 
sary node integers. 

l. Express the argument kas a binary number. 

Ze ion ate me array, shift the number y- £ places to the 


right and reverse the bits. 


To find out where the solid and dashed lines originate: 


1. Let the argument k be expressed as a binary number such 
tlarat pk = Rye! Se or Ko: Thus if k= 6, y= 3 then 
K, = 1, Kk = 1, Ky = (6) iche Kak Ko - 110), = 6) 509° 
Zs. Now in the ae array, node k has a solid line from Kk" 
node inthe £-1 array except that k ' is a one (1). 
= 
oe A dashed line originates from address k except that 
k is a zero (0). 
y-2 


To see how this works, apply the rules totwocases N ~- 4, N= 8 


and arrive at the tree graphs shown in Figures 3.1 and 3.2 respectively. 


| 
| 
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\ 
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\ <e 
| © 
ee ac gaComplex Addition 


ee Complex Multiplication 


N=4 yr2 


PIGURE 3. 1 


Tree Graph, N = 4 
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ee — .p»e-Complex Addition 


——______geee Complex Multiplication 


FIGURE 3.2 


Tree Graph, N = 8 
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Now verify the tree graph by looking at S(1l) = x, (1) in Figure 


Beal 
x,(1) =x, (2) + x, (3)W" 
x, (2) = x9(0) + x0 (2) Wo 
sate eal) x 9(3) Wo 3.5. 
Then, 


S(1) = x, (1) = X (0) - x 9(2)W* i (x9 (1) = x 9(3)W)W" 


S(1) = X (0) + x (1)W> + x 9(2)W* + x 9(3)W™, Bou 


which is in complete agreement with prior results. 


3.6 Advantage of the FFT 


The conventional method of performing Fourier analysis ona set 
of data points consists of summing a finite number of terms ina Fourier 
series to obtain each transformed data point. Each term added into the 
series requires a complex multiplication to obtain the term and a com- 
plex addition to include the term in the series. Thus, if N data points 
are to be transformed utilizing N terms of the series and one complex 
multiplication and one complex addition constitute one operation, then 
N¢ operations must be performed to accomplish the transformation. 
Likewise, a two-dimensional array of data points will require N 
operations. 

By comparison, the FFT requires N operations for each set 


of nodes and must have log,N sets of nodes which results in a total of 


1 








EZ 
Nlog,N operations for a one-dimensional transform. For the 
two-dimensional case AN‘ log,N operations are required. 
With the above results, the time savings factor (TSF) fora 


two-dimensional array of 64x 64 points is 


4N* log, N 4log,N j 
io = aoe S a = 70° Be Oe 
N N 


3.7 Problems in FFT Analysis 
Aliasing, leakage and the “picket fence’ effect are the most com- 
oT . 20 
mon errors encountered in utilization of the fast Fourier transform. 
The effects of these errors can usually be completely eradicated, or at 
the very least be reduced to acceptable levels, by insuring that the data 


is sampled at least twice during each cycle and by proper selection of 


data window. 


3.8 Space and Frequency Relationships 

PeuEouletmuEcicountened In performing either the Fourier trans- 
form or the inverse Fourier transform is the correlation of a sampling 
rate in the frequency domain to data point incrementation in the time or 
Special domain. In conventional Fourier analysis, the problem is 
resolved by the very nature of the series, in that the argument of the 
functions in each term of the Fourier series contains time and frequency 


explicitly as seen in equation 3. 8. 1. 


1 
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Cc 
x(t) = » om QJet nA ft Ke = 


n=-© 


A 
if 


The Fourier series does not demand the computation of x at any particu- 
lar t. The only requirement is that t lie within the data window 
-T/2 tS T/2 inthe asymmetric case and within 0 <t< TI inthe 
Syimimetric case when only N terms of the series are utilized. 

By contrast, time and frequency are nowhere to be found in the 
formulation of the fast Fourier transform. This can lead to confusion 
and uncertainty regarding the exact coordinates of a point in either domain 


unless careful attention is paid to equation 3.8.2. 


I on 
NAf NAw 





These relations are simple, but can mislead the user of the FFT. 

In order to clear up the question Surrounding this problem, con- 
sider a one-dimensional problem in which an array of data points are to 
be (1) transformed into the frequency domain; (2) multiplied by some 
other function; (3) then the product transformed back into the time domain. 


The problem may be stated as 


N-1 
S(k) = > x 9(n) eect N 3.8.3 


n=0 
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N-1 
j2nm nk/N 
(a) =o — STIG(k) eo / ance 
] N 
i=) 
with the following definitions (See Figure 3.3). 
T = width of data window 
N = number od data points 3.8. 


These independent choices dictate the following relations. 


a 
ore NE 
_ 2m (N-1) 
Aw = NT 
F = Spectral window = = 
a 
= —- oe 
Ax, Nol 3 


Then, assuming G is a continuous function, it must be evaluated 
in the interval -F/2=w = F/2 at N points spaced equally in increments 
of A w. This also assumes that G is a function of w rather than f, 
and this is quite often the case, e.g., the relative heating equations in the 
previous chapter which have been programmed utilizing the FFT and 
the results compared with a Gaussian quadrature transformation in the 


next chapter. 








CHAPTER 4 
UTILITZATION OF THE FFT ON RELATIVE HEATING PROBLEMS 
4.1 Introduction 


Diathermy, the therapeutic heating of deep body tissue, is an 
important part of physical medicine and rehabilitation. Research at the 
University of Washington in this field is concerned with the determina- 
tion of an optimum frequency and aperture configuration which will 
yield good heating in musculature without excess heating in the sub- 
cutaneous fat. This can be recognized ina relative heating pattern 
which shows low heating in the fat relative to the amount of heating 
observed in the muscle. 

The research basically has two phases. The first stage isa 
theoretical analysis of the problem utilizing plastic models of tissue 
and computer models of the electromagnetics. The second phase is the 
building and testing of diathermy applicators, once again using the 
already proven plastic and computer models. 

Chapters 2 and 3 have laid the groundwork for a computer model 
of the diathermy problem. The succeeding sections of this chapter 
will finish the construction of the computer model. Each phase of the 
research has a similar, but different computer model and each will be 


discussed separately. 
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4.2 The Theoretical Phase 

The first phase computer model for diathermy research consists 
of an executable computer algorithm which can evaluate the relative 
heating expressions of Chapter 2 for a theoretical aperture distribution, 
A suitable algorithm which can accomplish the transformations is the 
fast Fourier transform. 

Several versions of the fast Fourier transform can be found in 
the literature, some of which are written in Fortran and others 
are written in Algol, The algorithm of meewchi-: has been used in this 
research for two reasons; 

a.) its language is Fortran, 

b.) it can be implemented on the University of Washington 

computers without extensive revision. 

Some revision of Eguchi's program has been made to cut out the 
portions of his program which are not necessary in the diathermy prob- 
lem. Since extensive calculations must be done to build the arrays to 
be transformed and the amount of memory used is a determining factor 
in computer costs, features in Eguchi's program such as magnitude 
distributions, convolutions, gradients and the various output options 
are unusable and have been discarded since they obligate memory that 
can be more efficiently used. 


Figure 4,1 gives a flow chart for the use of the various sub- 


routines given in the Appendices and outlines the method of computer 


solutions necessary in the first phase of diathermy research. 





ot 


START 


READ Ax, AY, 
ZMAX, ZMIN, 
NCASES 


INITIALIZE 
AU, AV, UINIT, 


VINIT 





DO 82 
NMN=!,NCASES 


(yW READ FMC.EM, 
TM, EF TF ZI, 


> 
ow 


oH 


G ZERO 


SUMMING 


ARRAY 


FIGURE 4. 1 





Flow Chart of Theoretical Problem 








DO 8! 
INDEX =1,3 


V= VINIT 
DO 45 
M=1,64 


U=UINIT 
DO 44 


N=1,64 





FIGURE 4. 1 (continued) 





vV=Vv+Av 


PRINT PRINT 


lex 2 


radi}te: 


PRINT 
lez 


ADD INTO 
SUMMING 


ARRAY 


FIGURE 4. 1 (continued) 
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O-e 


PRINT 
RHP 


PUNCH 
FOR 
CALCOMP 


INCREMENT 


Zz 
oO 
! E 


YES 


Sel 


ibe 


FIGURE 4. 1 (continued) 
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4.3 The Empirical Phase 


Upon successful determination of an optimum frequency and 
aperture configuration, extensive tests must be made of any prototype 
diathermy applicators. 

Utilizing the methods of the first phase computer solution and 
having available a sufficient amount of empirical data from the field 
configuation at the surface of the plastic models, the fast Fourier 
transform can once again be used to determine the relative heating 
patterns throughout the volume of the planar model. 

A suitable method of solution for problems of an empirical 


nature is outlined in Figure 4. 2, 


4.4 Objective Restatement 


This research is primarily concerned with evaluating the relative 
merits of the fast Fourier transform to the solution of the theoretical 
relative heating problem versus the Gaussian quadrature method 
currently in use. Solutions for empirical problems will follow as a 
direct by-product of this research. The results of computations by 


both methods are compared in the next chapter. 
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> TART 


oF 
f 


= 
z 


READ AX. AZ, 
ZMAX ,ZMIN, 
Z\ 


INITIALIZE 
AU, AV,UINIT, 


VINIT,Z= ZI 





READ FMC.EM 
TER TRA,B 


FIGURE 4. 2 


Flow Chart of Empirical Problem 
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‘e ZERO 


SUMMING 
ARRAY 





DO 8s! 
INDEX=1,3 


qa V=VINIT 


DO 45 
M=1,64 


a U=UINIT 


DO 44 
N =1,64 
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FIGURE 4. 2 (continued) 
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PRINT 


EYP 





FIGURE 4. 2 (continued) 
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SUM INTO 
SUMMING 


ARRAY 
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PRINT 
RHP 


PUNCH 
FOR 


CALCOMP 


INCREMENT 


NO 


Oil 


YES 


STOP 


FIGURE 4. 2 (continued) 





Grey Par R™S5 
RESULTS 


5.1 Introduction 

Relative heating patterns (RHPS) as computed by the fast Fourier 
transform (FFT) and Gaussian quadrature (GQ) methods will be presented 
Mime@e succeeding sections of this chapter. An error criterion will be 
defined and the computed relative heating values will be compared ona 


point by point basis whenever possible, 


Pmeeeeenie Error Criterion 


Any digital integration method has an associated, inherent error. 
Quite often, the exact error cannot be ascertained; i.e., the bounds of the 
error are all that can be established. At other times, even the error 
bounds are difficult to establish. Inthe case of the FFT, Gentleman and 
pends” have proven the bounds of error due to round off, but nowhere is 
there a complete proof of error bounds for a fast Fourier transformation. 
No attempt will be made here to set error bounds. Instead, the FFT will 
be compared to the best available values from the Gaussian quadrature 
ime thied. 

Since both methods have inherent inaccuracies due to less than 
infinite integration, any error encountered in comparing results will also 
contain inaccuracies which are impossible to establish. For purposes of 


this research, error will be defined as the absolute difference in relative 
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heating values at the point (x,y,z) as computed by the Gaussian quadrature 
and fast Fourier transform methods. This error will be put ona percentage 
basis by dividing the magnitude of the error by the magnitude of the nor- 
malization point always has a value of unity in both methods, the point 
(0, 0, Zo) will always have zero error when comparing the two transforma- 


tion methods. This definition is then 
% ERROR = 100 (PGAUSS - PFFT) Bec. | 


where PGAUSS and PFFT are the relative heating values computed at 
(x,y,z) by the Gaussian quadrature and fast Fourier transformations 
respectively. 

To be an acceptable tool for computing relative heating patterns, 
the fast Fourier transform must compute values to within 5% ERROR, 


except at points which lie along lines with very steep slopes. 


5.3 Data Presentation 

Relative heating values and the associated error are tabulated in 
Tables 5.2 - 5.6 for several different frequencies. Each table is fol- 
lowed by CALCOMP plots which give a visual representation of the rela- 
tive heating patterns as computed by the FFT and Gaussian quadrature 
methods, im that order. Table 5.1 list the frequencies, associated 
dielectric constants, and loss tangents used in each of the five cases 
considered. East case considers a rectangular aperture (12 x l6cm.) 


in an infinite ground plane (See Figure 2.1). The relative heating 
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meeetas were Computed at Z = 0.1, 1.0, 2.0, 3.0, 4.0cm. with the 
fat muscle interface occurring in each case at 4, =e OC ii, 

Am briel explanation is in order concerning the CALCOMP plots 
which follow each tabulated comparison of the two methods. The first 
plot is a relative heating pattern as computed by the FFT. This is 
followed immediately by a transparency taken from the FFT results. 
This transparency overlays a relative heating pattern for the same case 
mescomputed by the GQ method. 

The transparencies unavoidably shrunk somewhat in processing, 
nevertheless they provide an excellent visual comparison of the two com- 
Butational methods. For example, examine Figures 5.3, 5.3a, 5.4 
which are the relative heating patterns for the 750 mHz comparison in 
iaoles>. a.) [he following agreements and discrepancies are readily 
identifiable in the overlay: 

a. There is very good agreement in the two plots for Z=1 

and 4 cm. 
Paiiiich7s=55 cra, curves difter by about 1% in the range 
3 cm. <x < 8 cm. 

c. About 2% difference can be seen inthe Z=2cm. curves 

in the range 1 cm. <x S$ 8 cm. 

deeeihe HET curve for Z2=0.1 cm. has some aliasing due to the 

truncated series, but shows good agreement with the GQ curve 
except at the edge of the aperture. 
Similar comparisons can be made for each frequency, but will not be made 


here. 





a, 


Table 5.1 
Input Data 
3 oe Ue é , Te 
400 MHZ Soe Ome 42 Sess. 0, 3745 
750 MHZ 51.54 , 5944 pee Blclens, 
918.8 MHZ bile . 488 Dat 5 faa 
foo MHZ. 49252 OSS op 76 - 1943 
2450 MHZ 47.42 Simo 059 Bee ao 0.1541 


Bees irequency 


dielectric constant in muscle 


m 
I 


m 

T = loss tangent in muscle 
m 

€ = dielectric constant in fat 

JT. = loss tangent in fat 
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Table 5.2 
400 MHZ 


Gomparison of _ RAPS as Computed by FFT and GQ 


—_— 


Z =O (tere = 0. 0 crm 4,2¢.0cm 
xX Pater | FP GAUSS % ERROR 
02:0 0.169388 0.169952 05056 
O.5 . 175147 SL Ose ws . 461 
0 . 171048 mee . 124 
I 5  lrooll . 175194 - 342 
ae O . 175963 pal ie LO ore 
Ze 5 . 185448 . 184532 O92 
aa 0 - 183975 . 190959 .698 
5.25 - 195594 . 198569 0.297 
4.0 - 194992 ero 51 eZ eo 
4.5 . 209297 S217296 0. 800 
50 . 209303 .228461 1.916 
5.5 ~ceiooal 241038 Poe 
6.0 . 228445 Pao oo6 2.714 
6.5 VZ5161556 won Osc Pee 0 
0 _ 209753 ogo 3905 
> Moo eoD | . 347174 1.432 
8.0 . 326463 - 437956 1.149 
Se: hom oo9 . 190027 3.844 
9.0 .090918 Pious no 3. 946 
9:5 "OGemt 5  LOoGs7 323.97 
O20 .051848 0.081288 2.944 
Zul. Or cune ye 0: 0 cm Z,=¢.0 cm 
0.0 . 091048 0.090634 O. 041 
0.5 . 091409 .091156 m2 5 
0) . 092435 .092712 n0ZcS 
ee 5 . 094167 . 095286 ey Z 
Za 7096532 . 098847 oo 
ZS . 099564 MUO SS 57 .379 
220 pL O Sako 5 -hOSt 71 959 
3545 . 107442 oLES0O9 oS 
4.0 Le 266 Pee 2 00 OF o3 
4.5 paelaic ¢ . 130148 1.242 
5a) we 37 53 . 138889 1.514 
be 5 08 7 2 . 148365 1.799 
be0) Opa 72co( 0.158386 Zi O 





ve 


nD .O OO =) = 6G 
Oo © Ol © Ure in 


— 


Dal 


Table 5.2 (continued) 


Owowowonrtnrt—rononn Pp hwWWN NR kK COO 


— 


WN NRK KK CO OO 


Ser @ Oi Se On © Ot ore Ol Om Ol em nt © Ol © ot © 


el Cy Orc) Cl 


foie wren  ¥-=. 0.0 cm 4, =4.0cm 
Peo P GAUSS % ERROR 
O. 143945 0. 168356 2.441 
. 148450 16ST 5 Zoe 
. 146659 . 178684 $2205 
S22 - L67 522 S589 
, LOZ lo5 . 143679 31051 
OG 3202 . 118478 32 526 
. 065198 097531 3.255 
0.052651 0. 080851 Ze620 
Z=92 0 cine 0.0 cm 4, = 2.0 cm 
1. OOO000 1. 000000 0. 000 
0.997064 0.996597 . O47 
. 988276 - 986425 ~ 185 
. 973692 . 969587 .411 
. 953406 - 946253 OF n> 
- 927548 . 916656 1.089 
. 896284 . 881081 eZ © 
-859815 - 839869 L995 
TOG © - 793416 2.496 
1 tec4l .742172 SOR 
s toler S | . 686660 32507 
.667224 . 627489 3.974 
. 609190 - 565400 4.379 
- 548269 50M 25 4.694 
- 485449 . 436549 4.890 
- 422356 . 372947 4.941 
. 361438 pols 0R 5 4.836 
. 305438 eat 4.585 
.256178 . 214073 4 oad 
.213941 . 176540 3.740 
tise 0.146014 57 
Zit= 3 einen ae 0. 0-Cm 4, =2.0cm 
2550851 0. 551602 0.125 
» 548625 . 549611 ~ 099 
~ 543457 . 543656 . 020 
. 534874 Soo o1o7 . 109 
a5c2 926 . 520097 eo 
. 507689 5027.10 498 
. 489268 0.481801 O. 747 
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Table 5.2 (continued) 
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Oo .O OO =) =) O Oo UF UT wd WW 


Oo oWonntnnrnurnh bHWWN DN RK KF OO OO 


own © Oro uUr1og oO Wi @ Ole oi 


So ou oOo Uo Uo UO Uo O1 OG urd oT oO Ur © 


Gea oe 0 ern. ¥ = 0. 0.eCm ZS a ae 

Perl P GAUSS % ERROR 

0.467811 0.457597 PeOZZ 
- 443515 - 430385 ie 35 
. 416643 . 400537 1-6 el 
. 387545 = Sievers) LY, 12903 
. 356670 - 334904 Cat 
. 324587 - 300397 2.419 
291979 . 265826 Ze outS 
- 259628 mc oel lO 2.752 
NOCOSS6 . 200174 2.816 
. 198829 meleos o2 2.800 
leialoe.1 . 144637 2.698 
. 146975 . 121814 2.516 
. 124916 . 102294 Zee be 

Of 05311 0. 085808 1. 950 

Jamar Oem Y = 20. 0-cm 2, - 2,0 cm 

0.300086 ORO TS) fee: 0.149 
7209063 . 300403 . 134 
. 296001 - 296889 . O89 
. 290926 O91 ORS SOs 
. 283882 Cog . 086 
.214937 . 272840 -210 
. 264192 . 260654 . 354 
Pe li oc - 246649 Pol 
.237889 oO US t . 683 
MOL tae -214170 Ossi 
. 206626 - 196335 1.029 
. 189871 SZ angoZ 
. 172844 . 159463 123136 
miso 92s . 141324 1. 460 
. 139465 . 123966 l-t5a0 
. 123766 BLOT? 56 L760: 
OD 0S 7 . 092960 1.608 
. 095389 Oman 21 eto 
. 082849 . 068066 1.478 
BO asi 057938 1. 345 

0. 060946 0.049218 dae 3 
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Table 5.3 
750 MHZ 


Comparison of RHPS as Computed by FFT and GQ 


— 


Zee Oe leectam. = 0. .0 cm B= 2.10 (eras 
x Pat P GAUSS % ERROR 
0.0 0.219999 O72 30121 LO eZ 
O.5 SEet 925 .-230939 0. 302 
0 IED OC t . 233344 DZ 
eS ee oD oD Soot biz . 174 
220 neo a © -241917 .814 
aS » 247462 - 247393 2007, 
3.0 .246062 pe oor U0 . 7104 
5. 5 ~coOULI2 . 258668 . 144 
4.0 220119 - 205705 .705 
4.5 .2710403 . 268169 woo 
5. 0 . 263880 eet le 79 . 800 
Se) ~-278808 21 5rl3 370 
6. 0 .269301 ~21ODa5 nL a, 
6.5 .291805 -283724 0.808 
toa . 280616 »o9391 5 iesS56 
ta - 345529 «320156 ZO 
8. 0 . 298479 SF 55C16 5.474 
8.5 .214655 . 144184 1. 047 
9.0 . 073904 07 E55 Oizo 5 
9.5 Ob2573 . 052983 0. 041 
0.0 0. 037884 0.041624 0. 374 
Zeal Orc mameye=.0. 0 Cm B= ae ee 
0.0 0.111816 O. 113543 0.173 
0.5 eae) . 114293 a heO 
10 ll SSO . 116457 bo 
LS . 118966 L977 35 SOeZ 
230 aClesioo > . 123896 . 026 
25 . 128689 . 128342 O35 
3.0 SS5Z28 . 13c635 . 088 
565 . Wor 741 . 136397 . 134 
4.0 . 140919 SOC! . 164 
4.5 . 142974 PL41132 . 184 
50 . 143856 . 141959 . 190 
ous) . 143801 . 141860 - 194 
6.0 . 142846 . 140888 . 196 
6.5 0. 140843 0. 138709 nel} 
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S22 00S OreSeen Oe Out Gor Oo Oo UO Our S 


Cie st ey Oi ey ti eS UT OS 


Zeon Oreim) Y= 0, 0 cm Z, = 2, 04a 
eee P GAUSS % ERROR 
0. 136485 0.133994 0.249 
.216978 1257603 oe 
.109112 a eOS sz . 380 
. 086657 . 083343 moo 
. 066680 . 065268 . 141) 
- 050719 - 051936 2422 
0. 037949 0.041896 Or 395 
fea ac Oemagey = 0.0 Cm 4, Orc ra 
1. 000000 1. 000000 0. 000 
0.993348 Oe IIS ToS . 043 
. 973627 OT D5Do a 
Oa S25 . 945266 ~ 374 
. 898147 . 904547 . 640 
. 844959 . 854473 02951 
mioon ll - 796589 1.288 
Ontos a4 ef aoe T Loe 
. 644953 . 664420 1.947 
2o7 594 - 593836 2.224 
. 498302 15022095 2.439 
. 426970 . 452710 Zor 4 
Zoo oy - 385434 Zao 3 
. 296831 - 322244 2.541] 
240923 . 264404 2.348 
. 192899 sucks Viel 2. O27 
- 153845 . 169734 1. 589 
» 124032 . 134680 15065 
mOzaZ 7 . 107445 O27 502 
087010 . 086560 OM05Z2 
OG ag od 0.070414 On 551 
Z=3 Cue) p= O..0 cm 4,2 2.0cm 
2320255 OF ao St O20Z5 
. 516642 75 looe835 . O01 
. 506000 7 506Go2 . 068 
. 488694 - 490472 Shs 
. 465339 - 468543 ao2.0 
. 436742 . 441603 . 486 
. 403861 . 410496 . 663 
P2501 to) Sonor Ow 0. 841 
. 329544 - 339601 1.006 
0.290378 sO beac eae 
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Table 5.3 (continued) 


(eee) ery = 70, 0. Cm aT 

ae eas P GAUSS 

0.251410 0.263964 
eO57 (4 . 226984 
. 178563 OTe 
. 146774 sl 5971S 
- 119219 erst liZ 
. 096394 - 06731 
Om 6552 . 086502 
. 064673 PO7 017 3 
. 054853 Ood Se 
. 047166 . 046785 
- 041558 0.038461 

Oram, | = 0.0 cm 4, 
me 70218 0.270042 
. 268288 . 268239 
. 262569 Baocood 
oo? oll wc 54162 
. 240690 . 242347 
NEC OLOS wee toce 
poor ao) .-211048 
. 188166 . 192560 
. 167699 rd 2 950 
. 146866 . 152844 
mlzo303 32655 
. 106861 ees 07 
. 088956 . 095864 
. 073124 . 079805 
. 059664 . 065816 
- 048675 BUSOU oo 
. 040046 . 0442 32 
. 033496 70S 6oN5 
. 028636 mua Oecd 
. 025039 . 024768 

0. 022295 0. 020567 


Seem Cm) 


= cm 


% ERROR 


me OD 
oc 
2330 
aio 
OS 
. 034 
7 OLS 
50 
6 Ott 
- 036 
0. 319 


OFM KF KF KF KE Ke 


0.018 
O05 
- 032 
. 090 
. 166 
oO oe 
- 347 
aoe 
Soa 
290 
Ooe 
66 5 
ooW 
. 668 
26 le 
5 om 
.419 
ZO 
eZ 
OZ 

Oe lis 
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Table 5.4 


OSs owl Z 
Comparison of RHPS as Computed by FFT and GQ 


Zi= 30) lserna  Y = 0.0 em a wc. 0 em 


p— 


x Eye ei P GAUSS % ERROR 
0.20 0.307665 0.319533 eee or 
O75 - 318697 . 320667 0.197 
1.0 . 312688 Oe OOo Leo 
ee POcColo | » 328736 0.055 
220 . 324581 . 334546 shel Tl 
2.5 » 341489 . 340525 . 096 
Ser) . 336246 . 346964 ORZ 
SD 251 G 5 1 . 350319 0.142 
4.0 . 342179 Ooo 1.114 
4.5 . 356449 35> 5055 0.136 
Su0) . 342905 . 356094 PacoshiG 
SD . 360407 oe Tt 0.319 
6.0 » 345085 mood ol 1.471 
6.5 Poor . 365824 0.976 
on . 359559 . 379465 ao 
Tao . 445967 - 413545 3.242 
8. 0 . 369174 -. 447123 Tee o> 
San . 147920 lo Z507 1. 541 
9.0 . 092290 080835 1. 146 
G5 OT O20 .059541 MAWo” 
OF 0 0. 054302 0.046853 0.745 
Zw lO era y= 0.0 Cm 4, 2c. 0 Cm 
QO. 0 0.144381 O. 146859 0.248 
O25 . 145761 . 148084 POOL 
ii} . 149585 wis b553 Pm we 
eed » 155296 . 156693 . 139 
ZO . 161878 . 162667 079 
225 . 168374 . 168539 m7. 
3.2.0 ~173745 . 173442 e030 
jan ol? (43 . 176752 . 066 
4.0 717901) . 178196 2OSZ 
4.5 Pitoe LZ . 177896 . 092 
50 . 177246 . 176304 . 094 
b> pelt I - 174003 ed 
6.0 . 172744 . 171344 . 140 
oe 5 . 169904 ots 25 es 
10 . 164432 61370 306 
> elioe 26> . 147447 . 482 
ee) 02130205 OF 22927 Ou ae 
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Table 5.4 (continued) 


Eee 


0. 105004 


] 


. 084526 
. 068273 
0.054134 


. 000000 


0.994130 


Otermaeey = 0. 0 ecm 2 


0 


- 976754 
5 yas, [3s 67) 
» 910485 
- 863800 
. 809766 
- 749643 
. 684639 
sole) 3 | 
» 544743 
- 472431 
. 400548 
yoga 70 
Goo > 
. 206994 
loans > 
. 119638 
709259 1 
. 074256 
OF 


062899 


7592 he 90 
elo 4 
- 508664 
- 492652 
sot ll35 
- 444872 
. 414672 
ooh ooo 
. 345638 


O2308315 


pOneminmetu= O20 Crm) Z 


PG ado 


0.095295 
Ue oo9 
. 058699 

0.047413 


l 


1. 000000 
0.994035 
. 976381 
. 947722 
moe 2U97 
. 861776 
. 807131 
. 746547 
. 681373 
. 612953 
. 542659 
. 471968 
402477 
. 335891 
1274031 
. 218845 
. 172206 
135132 
. 107008 
085977 
0.070018 


] 


0. 520860 
Poa 
5026 
491379 
- 469684 
. 443222 
- 412841 
. 379406 
oa og. 

0. 306731 


= 3. eine 


ere OC TY) 


= 2.0cm 


Poe ROR 


eae wil 
. 064 
Poot 
AG) (74 


oO Oo = 


0. 000 
OHO 
O57 
mS a 
5 UES: 
sauce 
. 264 
- 310 
aes 
eS 
iOS 
. 046 
aS 
eou0 
. 848 
OO 
. 442 
. 549 
. 462 
ali 
ye 


ORK | KH KO 


GnUS 
LOG 
Rolle 
ile @ 
» 145 
OS 
oes 
ae a Je) 
. 188 

02158 
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@ 30 0 OO =) 15" Oo OT oi 


DOO WONNTAA UNA HWBWN NH kK CO OC 


ae) OC Ole oO ea at aa © 


© 


Se oS OVO wm Oo Clo oO Oro U1 co UI © oO 


Zao eCtmiane =O, 0 cm Bo 


Lea 


0.270145 
goo l 93 
. 194862 
MEDI e SS 
Mico Leo 
. 100900 
Ome oc 9 
OG 12 17 
. 048567 
US 9659 
. 034219 


Ocm Y = 0.0 cm Z 


0. 
men 0639 
Zone 
cool | 
. 244061 
ea} 
polcso3 
aoe 
. 174496 
- 154269 
Loo aS 
~ 114105 
~ 095275 
PUR oUS4 
SUG2 ocr 
. 049886 
nO soc 
sOslis5 
BOZO 5s 
7020963 


ae 


pic 516 


OeUlec 43 
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Table 5.4 (continued) 


PG AuUSS 


0.269174 
HoS3cUU0S 
. 196241 
. 162945 
als oleic 
5 IG See 
. 086485 
069107 
7056560 
. 046257 

0. 038083 


On2 (15 96 
-coom ls 
- 264150 
CODD E42 
2243057 
. 228344 
ee LOO, 
lS 0ne 
le 
miss 5a. 
TiS 5655 
- 114382 
~O95325 
Un ono 
pUCSO0% 
RUD oO56 
ot Oil 
. 035806 
02922 
. 024413 

0. 020317 


= 2 Ose 


=O Cm 


% ERROR 


0. 097 
OT 
cl 3S 
=oUG 
0 
. 664 
2196 
. 849 
Stee 
. 640 

0. 386 


0, 092 
ODS 
. 094 
One 
> LOG 
ol0O4 
. 106 
‘OS 
O02 
2069 
- 030 
, 028 
. 103 
> EZ 
cise 
tS 
439 
paOe 
- 432 
- 345 
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Table 5.5 
1500 MHZ 
Comparison of RHPS as Computed by FFT and GQ 


ua) items = 0). 0s Cm 4, =42.0cm 
xX er P GAUSS % ERROR 
0. 0 0.769940 0.791965 ee Oe 
On 5 - 797484 on eZ 0.238 
ie 0 Siac ol Sous 2 (op ON Sye, 
eS mole ot Se 4 Shey | 0.556 
ZO . 800998 .819285 1.829 
2. 829700 . 821610 0. 809 
5.0 -sv05l2 Be VA 1. 941 
ono . 824841 . 817286 0.756 
4.0 at oo c20C , Olt aot 2.435 
4.5 oo O36 ~Og5a01 0. 688 
ore .805427 . 834203 ZO to 
S25 . 860066 847315 ls Cre 
On) . 827499 . 857974 3.048 
6.5 -892359 . 862831 21 993 
a7) - 821465 - 863583 4.212 
> - 957599 . 876068 S, 158 
Caw) . 642953 . 781976 3.902 
SD . 164369 . 155495 O. 887 
9. 0 0% 7 bo3 . 083483 OFGec 
9.5 . 050834 . 061547 1.071 
10. 0 0. 044408 0. 050588 0.618 
em oli Cmy erat OO Cimaa Bo = 2. (chee 

Oe 0.371453 0.375005 07 355 
Om5 . 376461 - 379559 oe 
ea) . 389622 ole eed 
hed . 407001 . 407720 Ors 
Zi2 0 . 423560 »423174 . 039 
Za 4361 71 , 435057 - Le 
300) . 443940 . 442821 Sle 
S745 - 449051 . 448231 O32 
Ae . 454231 . 453964 FOZ, 
4.5 »461729 . 461654 . 008 
Ba) . 470444 . 470326 = Oe 
ba . 476507 - 475855 , 065 
oe) . 472567 471472 aad) 
b> » 450246 . 448732 bol 
eee . 400007 2596953 . 145 
7a . 316297 2315] 45 LS 
8. 0 0.212149 Omecileeen stan 0. 063 
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ables se. > (continued) 


lOve = Os 0 cm 4, =¢.0cm 


hee P GAUSS % ERROR 
OF PZ28503 Opts 34i75 0.497 

. 081896 Ora O50 

. 059079 . 069619 1.054 
0. 049324 0.056038 Oho val 


ZO ser or = 0. 0 em Z,=¢.0cm 


1. 000000 1. 000000 0. 000 
Zo ay lL -U2Z2908 . 024 
1.090709 1. 089751 096 
it 96503 get me) 
io 23595 1. 325446 , 344 
1.470068 1.465576 . 449 
77 090 io 92>2° . 454 
1634959 leoo2ues - 284 
Pear il So) 1.712660 rt 
loge a9 lor Ulios Oran 
13550 le 5o2269 1. 466 
1.347428 1. 369424 2-00) 
iat Lola 0 Ca oo 
O25872545 O07 90lG37 ZO 
. 646177 or aOne Zev) 
- 458807 - 476519 Poa? | 
= 321036 . 327093 0. 606 
ye D0 5 / .224164 O33 9 
eS 53555 1.424 
- 133854 90 1. 606 
0.102661 On Owe lao: 


Ze 3 Overman eo 0. 0Fen) 4a 0 Cm 


mm HW WN DYN — — CO OO 
whee UT CSV Oe oy 2 U1 © 


] 


0.445169 0. 443357 0.181 
. 454701 . 452869 mos 
- 482492 . 480589 5 SN, 
ToCooUt Soca ome . 204 
5.9706 ot 1 lew roe 
BOSiOo oe . 634624 S225 
. 686838 . 684962 . 188 
5 NO rete i lhemec . 086 
miceo OO - 725647 . 100 


0.697849 O02 701516 02308 
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Table 5.5 (continued) 


Zo one -Y = 0,,0-cm 4, 
Poe P GAUSS 
0. 639040 0.645835 
noo oS - 564262 
- 455316 . 467056 
~ 354539 . 366590 
. 263993 . 274248 
- 191319 . 197768 
. 138623 mio! 
O31 30 . 099480 
Oro 4 oo . 072423 
- 062145 . 054623 
O. 047676 0. 042668 

feo eer YY = 0.0 Cm 4, 

O. 198890 0-197 535 
meULI 2S .201401 
-213910 m2 26 On 
,231335 Pe oUlS! 
PZOCOoT Peo 51664 
275344 274490 
2.94464 - 293935 
, SUB OS . 305781 
1 05505 306291 
BMagisoS - 293457 
. 264338 SConelS 
ReZ21 505 1252026 
. 185905 . 191090 
. 144990 . 150078 
. 109152 mel 3250 
. 080784 . 083091 
. 060134 2060205 
. 045859 » 043822 
. 035893 #032519 
BUZSZ17 . 024815 

0.021741 0.019488 


=e 0 4erMm 


= 2205emM 


% ERROR 


0 
0 
} 


om 


0 


0 


-679 
or 
.174 
~ 205 
025 
. 645 
Se 
_ 3005 
7 t02e 
ao 
DOt 


= le 
ol 
Peles 
1 ree 
108 
. 085 
O55 
. 004 
095 
. 209 
: 20 
» 452 
Soe 
09 
. 410 
coil 
. 007 
. 204 
oo 
. 346 
ee. 
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Pale 5.76 
2450 MHZ 
Comparison ot RHPS as Computed by FFT and GO 


Ze— aici © = 10.0 °em 4, 24.0 cm 


— 


x Pee tT P GAUSS % ERROR 
0, O 0.366811 0.377543 OS 
75 . 378686 . 378658 0. 002 
1.0 - 370445 . 380583 1.014 
i> . 381763 . 380782 0.098 
2.0 - 368881 . 378765 0.988 
Zz 5 SOLO L » 300030 0.145 
5a . 366957 ~ 377359 1. 040 
5.5 ~o6c005 . 380061 0.194 
4.0 . 370698 - 382443 1.174 
4.5 . 384621 ~SOCOrc 0.195 
a () . 367049 1901338 1.429 
55 oes cl - 380048 07 5108 
6. 0 . 363116 7 S107 10 1 56D 
6. 5 - 386408 - 3715910 1.050 
tO . 351085 eo 7co5 | 2.105 
ie) . 407494 25 (O50 Se Ge 
8. 0 oct | . 335005 6.433 
8.5 - 056666 ORG CL 0. 444 
9.0 015730 - 015568 0.016 
9.5 . 004814 . 006221 0.141 
0. 0 0. 004538 0. 004280 0. 026 
Zea ole Onetimet e—. ). 0 °C im 2, =¢-0 cm 
0.0 0.468521 0.471441 0.292 
0.5 . 468444 4.71251 meow 
0 - 466835 - 469570 oS 
cD - 461953 . 464718 Hoo 
220 . 454289 - 457362 oO 
Ze, 5 . 447671 - 450967 930 
3.0 - 445632 . 449004 od 
3:5 - 448780 - 451928 wok 
4.0 - 454143 457113 moot. 
4.5 - 458535 . 461291 cao 
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5 
Saya Analysis of Results 

The different ranges of % ERROR and the number of points lying 
within the various ranges are tabulated in Table 5.7, Examination of 
the six compared points which lie outside the range of acceptability 
reveals that each of the points has three factors in common 

a. all occur near the edge of the aperture, 

Peeealioccur when Z = 0.] cm, i,é€.; Zis small, 

c, all occur inan area where the relative heating pattern 

has a steep slope. 

It was mentioned in section 2,6 that points in the area defined 
by (a) and (b) above may need corrections to account for less than infinite 
integration, These corrections have not been made in the FFT, 

Since only 1.14% of the points compared fail to meet the 5% 
ERROR criterion and since all of these points are in an area of steep 


slope, the FFT is, therefore, an acceptable computational tool. 


Tee io 1 


Ranges of % ERROR 


% ERROR 
Greater But Less Number Percentage 
Than Than of Points in Range 
10% -- 2 0. 38 
5 10% 4 0,76 
* 5 8 Toe 
3 4 17 3,24 
2 3 34 6.48 
1 2 71 15,05 
0 1 381 logol 





Total 525 100,00% 
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Se 5 Cost Analysis 


In Chapter 3 it was stated that the fast Fourier transform is an ex- 
tremely efficient method of performing discrete Fourier transforms, In 
this light, two further points are noteworthy in a comparison of the FFT and 
the Guassian guadrature, 

First, all of the computed results in this chapter were performed on 
the CDC 6400 using the RUN (S) compiler, Onthe basis of central processor 
time only, an average time was taken, over several thousand points, for the 
compution of one point with the following results for the two methods, 

Gaussian guadrature = 1.9306 sec. /point Sol 

Fast Fourier transform = 0.01499 sec. /point. Bou 


This yields a time savings factor (TSF) of 


0.01499 1 
1.93062 }128.8 





TSF = 


which is 76% of the theoretical TSF derived in Chapter 3, 

The second point which makes the FFT not only acceptable but also 
desirable is that the FFT computes 4096 points for each x-y plane investigated, 
These points are, for a given Z, in the range 

-l6cm <x <1l6cm -l6cm <y <1l6cm, 

By contrast, the Gaussian guadrature computes only 21 points, for a given 
~, inthe range 


yas 0 <x <l0cm 





a 


advantages of the FFT are, therefore: 
reasonably comparable accuracy, 
greater time savings, 

mene «eonplete data, 


no restrictions on aperture distributions, 
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GCONGEUSION 


6.1 The FFT in Current Applications 

Current development of a safe, efficient, direct contact diathermy 
applicator for use in treatment of musculoskeletal disorders assumes, 
for analysis of the electromagnetic fields, a linear superposition of an 
infinite set of plane wave solutions, The general solutions of the scalar 
wave equation thus obtained develop a complete set of Fourier transform 
pairs which completely specify the electromagnetic field components 
throughout the volume of the human tissue model. 

In the past, the Fourier transforms have been computed by the 
Gaussian quadrature method. Economical use of the Gaussian quadrature 
imposed definite restrictions on the solutions thus obtained: 

a.) All theoretical and empirical solutions were forced to have 
complete symmetry about the x and y axes thus enabling 
integration in only one quadrant of the u - v plane. 

b.) The economics of computer time dictated that the relative 
heating be computed only along the positive y-axis in any 

—~ desired z-plane. 

c,) Extensive corrections were required near the edges of 
the aperture. 

These restrictions obviously limited the generality and completeness 


of any relative heating patterns thus obtained. 





lies 

The fast Fourier transform package developed in this research 
completely circumvents these severe restrictions. Integration is per- 
formed over all four quadrants of the u - v plane, thus eliminating the 
need for symmetry. In addition, the efficiency of the fast Fourier trans - 
form allows the computation of 4096 points in each z-plane rather than 
just a few points along an axis of symmetry. Further, even though no 
corrections are included in the fast Fourier transform package, nearly 
ninety-nine percent of all points calculated agree, to a reasonable 
degree, with answers from the Gaussian quadrature method. The FFT, 
on a point for point basis, accomplishes all of this in two orders of 
magnitude less time. Two minor disadvantages are encountered in the 
use of the fast Fourier transform. Some aliasing is encountered in 
z-planes near the aperture and somewhat more actual computer time 
is required for the computation of each z-plane. Nevertheless, these 
disadvantages are overwhelmingly outweighed by the completeness of 


the relative heating patterns obtained by the FFT. 


6,2 Future Application of the FFT 
Although the FFT has at this point only been adapted to the 
analysis of relative heating in planar models of human tissue, the 
. ZS 
package developed in this research should be extremely useful to Ho 
in analyzing relative heating in cylindrical models. 


Quantitative studies in the adverse effects of electromagnetic 


radiation will result in the establishment of radiation safety standards. 
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These quantitative studies will require an analysis of the electromagnetic 
fields in both planar and cylindrical models of human tissue. The fast 
Fourier transform package developed in this research can be easily 
adapted to this problem since no symmetry is required in the arrays to 
be transformed, 

These and other problems requiring Fourier analysis in two 
dimensions can be implemented by replacing the SCRPEX and JAMXX sub- 
routines by appropriate subroutines which build the two-dimensional 


arrays to be transformed. 





APPENDIX A 
SUBROUTINE NUMGEN 


This subroutine utilizes that portion of meveniuess program which 
generates the weighing factors, wos for each node in the transform 
and calculates the size and beginning and ending values for the DO loops 
used to transform each array. The arrays built by NUMGEN are called 
byt AST in Appendix B. 

If the number of data points is consistent for all uses of this 
package, NUMGEN need only be used the first time the package is used. 
In the first run, the arrays built by NUMGEN can be punched out and 
read in as data in subsequent runs. This will result in slightly less 
central processor time and somewhat more peripheral processor time 


which is degirable on the CDC 6400. Data transmission is via labeled 


common, 
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A Pee NbpIX B 


SUBROUTINE FAST 
: ae : ee 
Subroutine FAST utilizes that portion of Eguchi's program 
which calls the subroutine FFT from the bit reversed address arrange- 
ment, FFT is the actual workhorse that performs the Fourier trans - 


formation, Data transmission is via labeled common. 
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SUBROUTINE SCR PEX 


Subroutine SCRPEX calculates the Fourier transform ofa 
theoretical aperture distribution for a givenu and v. In Chapter 2, 
the final set of equations for the relative heating problem contained a 
¢ which was the Fourier transform of the aperture distribution at 
z=0Q . This subroutine is entered into the package as a separate 
entity to enable a potential user to change aperture distributions without 
rewriting the JAMXX subroutine. Data transmission is via labeled 


Coram on. 
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me Phe. Dp 
SUBROUTINE JAMXX 


JAMXX computes the appropriate value in the u - v plane for the 
array which is to be transformed. JAMXX requires an input from 
SCRPEX which, for a given u and v, corresponds to the Fourier trans - 

: it 
form of the aperture distribution. This subroutine by McDougall 
incorporates all the theory of Chapter 2 except the corrections and 
&, the aperture distribution transform. Data transmission is accom- 


re] 


plished via labeled common. 
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ore END x 


SAMPLE PROGRAM 


The sample program was used to calculate the relative heating 
patterns in Figure E.1. It utilizes all of the subroutines in the previous 
appendices in the manner outlined in Figure 4.1. Figure E.l isa 
CALCOMP plot of some of the wealth of information in each transformed 
array. This figure shows the relative heating pattern for various y's 
at the fat-muscle interface for 2450 mHz, It should be mentioned that 
the Gaussian quadrature, for the same problem, obtained only the 
relative heating pattern for y=0, x >0. Data transmission with the 


subroutines is via labeled common. 
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